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We extend a tight-binding total energy method to include /-electrons, and apply it to the study 
of the structural and elastic properties of a range of elements from Be to U. We find that the 
tight-binding parameters are as accurate and transferable for /-electron systems as they are for 
d-electron systems. In both cases we have found it essential to take great care in constraining the 
fitting procedure by using a block-diagonalization procedure, which we describe in detail. 
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I. INTRODUCTION 

In order to use atomistic modeling to predict ma- 
terials properties it is necessary to accurately deter- 
mine the forces between the atoms in a solid. Most 
of the established molecular dynamics simulations in- 
volve either pair potentials, embedded-atom potentials, 
or, for covalent materials like Si, some classical poten- 
tials with additional bond-angle information. Unfortu- 
nately, these methods have often proved inadequate for 
complex transition-metal and /-electron systems, which 
have important bond-bending forces that seem difficult 
to capture in any simple way. A promising approach for 
such systems is tight-binding (TB), since it automatically 
builds in the quantum-mechanical bonding that conven- 
tional first-principles local-density approximation (LDA) 
or gradient corrected (GGA) band-structure calculations 
can very accurately determine. 

In this work we report on our progress to extend the 
recent TB total energy model that was developed at the 
U.S. Naval Research Laboratory (NRL); this has been 
successfully applied tc^stmijconductors and both simple 
and transition metalarBaClij. We have extended and 
tested the NRL-TB scheme to include /-electrons, which 
makes possible the exploration of materials with impor- 
tant / characteristics, such as the light actinides. We 
also describe our own specific methodology for our fit- 
ting procedures, which is based on the NRL procedures, 
that we believe should produce parameterizations that 
are highly accurate and should be well suited for molec- 
ular dynamics simulations, although we report no such 
simulations in this paper. In the rest of the paper, we 
first briefly recapitulate our modifications to the NRL 
approach, and describe our experience in obtaining high 
quality TB parameters. We test this approach for an s-p 
bonded material (Be), c£-electron materials with nearly 
filled and partially filled d-bands (Cu and Nb), and fi- 
nally for a partially filled /-band system (U) , an element 
of significant technological importance, which has com- 
plex structural properties. 



One of the major obstacles in atomistic methods for 
materials modeling is the question of transferability, 
which arises from the semi-empirical nature of the atom- 
istic forces, which involve fitting parameters to either ex- 
perimental data or to theoretical calculations. The prob- 
lem of transferability is the question of how well such 
fitting parameters will work for geometries of atoms that 
are different from those to which they were fit. This 
question is especially important for molecular dynamics 
simulations, where the atoms are free to move wherever 
the atomic forces push them. It is also needed for de- 
termining ground-state structures — the structure that 
will minimize the total energy of the system — or for cal- 
culating various defect structure (vacancies, interstitials, 
and more complex defects) that are expensive for first- 
principles methods (LDA calculations). In this paper we 
will show that the TB method works for a wide variety of 
different atomic environments (e.g., for crystal structures 
with very different numbers of near neighbors) and over 
a large atomic volume range. 

Another qualitative feature of atomistic potentials that 
is required for accurate molecular dynamics and other 
atomistic simulations is to parameterize self-consistency 
or local environmental effects. To understand what 
we mean by this statement, consider first-principles 
electronic-structure calculations for the total energy of 
a solid as a function of the atomic positions. Such cal- 
culations depend on determining a self-consistent poten- 
tial for the electrons for every atomic structure consid- 
ered. An electronic potential for any given atom depends 
importantly on the local arrangement of nearby atoms. 
This implies that any force parameterization that suc- 
cessfully mimics the electronic-structure calculation also 
senses changes in the local atomic environments. This 
concept has been a major driving force for adding lo- 
cal environmentally sensitive terms in atomic force mod- 
els. The NRL-TB parameterization also contains such 
effects for onsite terms (see below). One strong test 
of local atomic rearrangements is the number of near- 
neighbors, which should strongly perturb the local charge 
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density around any given atom and hence exacerbate self- 
consistency effects. For this reason we have focused our 
fitting procedures for TB parameters on crystal struc- 
tures with widely varying numbers of nearest neighbors. 
As we will describe below, the NRL parameterization ap- 
pears to have adequate flexibility to successfully pass this 
fairly stringent test. 

All empirical force models use some procedure to fit 
their parameters (usually some variation on a least- 
squares fit to either experimental or theoretical data) . As 
more parameters arc added, perhaps to better describe 
bond-bending forces or approximate environmental ef- 
fects, overall properties such as the total energy often be- 
come less sensitive to the parameters. What is actually 
worse is that different parameters can sometimes com- 
pensate for each other so that parameter space becomes 
a vast multidimensional landscape where the minimizing 
function often has large numbers of weak minima and 
it is difficult to find the physically relevant true global 
minimum that provides maximum transferability. Since 
there is often a choice of different parameterizations, it 
can be difficult to know whether the lack of success in 
fitting properties is due to a parameterization that has 
inadequate flexibility or a failure to find and lock onto 
the best choice of parameters when wandering through 
the enormous parameterization space during the mini- 
mization process. 



In our TB fits, we have found it to be enormously im- 
portant to add constraints to the fitting procedure in or- 
der to force the parameters to be near the physically rel- 
evant portion of parameter space. Once this is done, any 
minimization procedure seems to improve the minimizing 
function as well as the transferability. Using physically 
motivated constraints appears to remove the huge degen- 
eracy of multiply compensating parameters and to help 
focus the parameterization onto the desired solution. In 
the case of TB, as we will discuss below, the key con- 
straint is symmetry. If used properly, it appears to force 
parameter space into highly transferable parameteriza- 
tions. 



II. TB METHOD 

Theoretical justification for the tight binding method 
rests upon the division of the total energy of the system 
into a repulsive pairwise term and contributions from the 
valence band structure, 

Et — E rep + Eb, (1) 

itself an approximation to the Harris-FoulkesiS func- 
tional, which caiiipbe derived from Kohn-Sham density 
functional theoryffl (DFT), 



S[n] 



^2 _ E H[n v ] - J n v (e xc [n v ] - V xc [n v }) + ^^^Z va Z vb / |R a - R f 



(2) 



a b^a 



In Eq. g, the subscript v denotes the fact that we are 
dealing only with the valence electrons whose density is 
given by n„, and the first term on the right hand side 
is the band term - a sum over eigenvalues arising from 
a Schrodinger-like equation. Eh denotes the Hartree en- 
ergy, and e xc and V xc the exchange and correlation energy 
and potential, respectively. A simple decomposition of 
this functional into pairwise and bonding terms (Eq. [I]) 
can be justifiednQ when terms involving more than two 
atomic centers are ignored (see Eq. ||] below), and the, 
input charge density is atomic-like. The Slater-KostertU 
method then consists of solving the secular equation, 
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for the single particle eigenvalues and orbitals, under the 
restrictions: terms involving more than two centers are 
ignored, terms where the orbitals are on the same atomic 
site are taken as constants, and the resulting reduced 
set of matrix elements are treated as variable parame- 
ters. Note that Eq. ^ includes the overlap matrix, S, to 
take into account that the basis functions need not be or- 
thogonal. To better illustrate these restrictions, we write 



the Hamiltonian including the labels for orbitals having 
generic quantum numbers a, (3 localized on atoms 
where the effective potential is assumed to be spherical, 
and can be represented as a sum over atomic centers, 



H, 



^ 2 + zZ v f 



0,3), 



(4) 



which wc further decompose into "on-site" and "inter- 
site" terms, 



Hai,/3j = e a S a p5ij + E ai ^j(l — Sij), 



(5) 
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where the on-site terms, e a , represent terms in which two 
orbitals share the same atomic site, and, for j =/= i, 



n 

x J drip a (r - R„ - bj) Hipp (r - b,-) , (6) 

^ = E elk,(R " +b ^ b,) 

n 

x drip a (r-Rn- hi) ipp (r - bj) , (7) 



are the remaining energy integrals involving orbitals lo- 
cated on different atomic sites. In Eqs. we have 
used translational invariance to reduce the number of 
sums over bravais lattice points {R n }, and the bj denote 
atomic basis vectors within the repeated lattice cells. Or- 
thogonal TB (OTB) treatments do not use Eq. 0, and 
thus have fewer parameters, but we have found a non- 
orthogonal (NOTB) approach ±q be more accurate, con- 
sistent with previous TB worklld. All of the results that 
we report in this paper use a NOTB model. Note that 
terms which have both orbitals located on the same site, 
but the potential on other sites have been ignored - these 
contributions are typically taken to be "environmental" 
corrections to the on-site terms, and are not accounted 
for in the usual Slater-Koster formalism (although such 
corrections have been explored, for example, see Ref. 11 
and references therein). The two center approximation 
then consists of ignoring additional terms in the inter-site 
contributions in which the effective potential does, not lie 
on one of the atomic sites. This approximation!!] is not 
necessarily very accurate (see, for example, the compar- 
ison between two and three-center fits in Ref. but it 
is quite often used due to the enormous simplification 
of the overall TB method; perhaps one can view TB as 
having the right functional form and the choice of pa- 
rameters allow some correction for any errors due to ne- 
glected terms. Once this approximation has been made, 
well justified or not, the inter-atomic (i ^ j) matrix ele- 
ments reduce to a simple sum over angular functions and 
functions which depend only upon the magnitude of the 
distances between atoms, 



H a i,pi = hui m (rij)Gii' m (Cli t j), (8) 
S a i,pj = y^su> m (rij)Gw m (flij), (9) 



where we have now adopted the usual convention of us- 
ing the familiar l,m angular momentum quantum num- 
bers, and the axis connecting the atoms is the quantiza- 
tion axis. The basis set used for the a and (3 quantum 
states are the cubic harmonicsE3 whose functional forms 



are given by (with appropriate normalization factors) 
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where /;(r) = l/r l , and |±) denotes the spin-state. The 
selection of these particular functions is not accidental, 
as they are chosen to specifically possess the various ir- 
reducible representations of the cubic point group Oh- 
In principle we are free to choose any set of orthogonal 
functions as our basis, but it is beneficial to choose the 
set that best reflects the symmetry pmperties of the sys- 
tem under study. Other choices existtj, but they do not 
provide nearly as transparent a framework for modeling, 
fitting and understanding the various parameters. 

The Slater-Koster tables for the st£d 5 matrix elements 
can be found in standard references^, and, we have used 
the tabulated results of Takegahara et al£3 for the addi- 
tional matrix elements involving /-electrons. Typical TB 
applications are then reduced to using TB as an interpo- 
lation scheme; the matrix elements {hu'm, Hl'm (if used) 
and e a ) are determined by fitting to ab-initio calculated 
quantities such as the total energy and band energies. 
The NRL-TB total energy method was an innovation in 
that the energy bands used in the TB fits were shifted 
such that the integrated band energy was the total en- 
ergy, 

e' b = {e b -E b + E tot )/N v , (10) 

where £& are the unshifted energy band values, E b = 
Y] b e b the total band energy, and E to t the computed to- 
tal energy. The NRL-TB method imposes a simple func- 
tional form on the inter-site matrix elements, 

hw m {r) = (a Wm + b Wm r) exp (-c 2 ,, m r) f c (r), (11) 
swm(r) = (awm + b Wm r) exp (-c 2 ,, m r) / c (r), (12) 

where f c = 1/(1 + cxp(2 * (r — r ))) is a multiplicative 
factor included to ensure a smooth cutoff with increasing 
distance. Our applications have used ro = 13.5 bohr 
radii for Cu and Nb, tq = 10 for Be, and ro = 12.0 for U. 
The on-site terms include the environment dependence 
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necessary to allow the parameterization to be applied to 
other structures not included in the TB database, 



(13) 



The p function is intended to be a crude measure of the 
atomic density (providing a correlation to the number of 
atomic neighbors), 



P = X] cxp (~^ r ij) fc( r ij)> 



(14) 



where is the interatomic distance. Note that, for a 
NOTB model in our sp 3 d 5 f 7 basis, this TB model has 
1 + 16 + 2 * (3 * 20) = 137 free parameters, whose deter- 
mination is the key to an accurate representation. 



A. Fitting the Parameters 

In order to determine such, a.ria*ce set of parameters, 
the usual NRL-TB approacrJlffloO involves performing 
a non-linear least squares minimization, fitting to the en- 
ergy bands and total energies calculated from density 
functional theory. The energy bands are fitted over a 
set of points in the irreducible wedge of the first Bril- 
louin zone (IBZ). In practice we have found it absolutely 
necessary to take maximum advantage of the symme- 
try present in order to impose constraints on the fitting 
process. We exploit the symmetry of all possible high- 
symmetry points and lines in the first Brillouin zone in 
order to reduce the possibility that the fitting process can 
easily mistake the ordering of the energy bands. We have 
decomposed the TB wavefunction at each high symme- 
try point and line in terms of the symmetry-adapted TB 
basis functionsEj, which allows us to block-diagonalizc 
the Hamiltonian and overlap matrices, and determine the 
eigenvalues corresponding to each irreducible representa- 
tion, which avoids possible confusion as to the band or- 
dering. This process is illustrated schematically by Fig- 
ure [|. Figure [l] shows the perils of fitting using a simple 




X 



FIG. 1: A schematic representation of the band ordering 
problem when trying to obtain high quality TB fits. Here the 
actual bands correspond to two different irreducible represen- 
tations, labeled Ai and A2. The TB fits, when the symmetry 
is ignored and only the energy ordering is used, can miss es- 
sential features, like the band crossing pictured here. In this 
figure the A direction in the Brillouin zone is the abscissa, 
while the ordinate is increasing energy. 



TABLE I: Basis functions used in the block-diagonalization of 
high-symmetry point X in the fee lattice. Only one of the two 
partners is shown for the two- dimensional representations. 



Xi 
X 2 
X2' 

x 3 
x 3 - 
x 4 , 

X, 

x 5 < 



Aig 

Big 

B 2g 

B2u 

Bzu 
E a 



{\s),\ds)}, 
{\d4)}, 
{|/i>}, 
{Mi)}, 

i>3),|/4)}, 

{(|pa> + b3»A/2}, 

{|Pl>,l/2),l/5», 



energy ordering scheme for the bands - unless the irre- 
ducible representation is identified and constrained for 
each band, the actual fitting errors are misidentificd. 

Let us describe the scheme mathematically and follow 
up with a simple example. Along symmetry lines and 
planes in the Brillouin zone, there exist symmetry oper- 
ations for each k that form a group. For each such group 
the irreducible representations T p can be found, and can 
be used to block-diagonalize the secular equation, 



D{p) Jl = (<jf>\H\^)-E(^) 



(15) 



where D(p) now takes on the "block- form" , and the sym- 
metry reduced eigenvalues can be found by solving 



D(l) 






D{2) 
D{3) 



(16) 



or simply |D(p)| = for every p. 

The resulting set of symmetry adapted basis functions 
are specific linear combinations of the original cubic har- 
monic basis functions (Eq. [n]), which can be found, in 
practical terms, by applying projection operators of the 
irreducible representations to the original basis functions. 
In our fits, we require that only the correct symmetry 
adapted tight-binding basis functions can be used for fit- 
ting each energy eigenvalue for high symmetry k points 
and directions in the Brillouin zones for the various cubic 
crystal structures. In this way we fit not only the en- 
ergy eigenvalues, but also we restrict the eigenfunctions 
as well. One must then also provide the TB fit with 
a database in which the energy bands are broken down 
by their irreducible representation at each high symme- 
try point and direction. In practice this means that the 
linear combinations of TB basis functions that have the 
irreducible representations at each k must be determined. 
These linear combinations then "block-diagonalize" the 
Hamiltonian and overlap matrices as described above. 

The use of symmetry in constraining the fitting process 
has previously been developed by PapaconstantopoulosQ 
for sp 3 d 5 basis. We have found our procedure for us- 
ing symmetry constraints to be essential in limiting the 
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parametric phase space and thereby obtaining transfer- 
able TB parameters. Here we extend the scheme to in- 
clude /-electrons. We wish to strongly emphasize how 
important this procedure has been in determining an ac- 
curate fit, especially with the large parameter set used 
when fitting to an sp 3 d 5 f 7 basis. As an example, in the 
fee lattice structure, each energy band at the point X 
[with cartesian coordinates (0, 2ir/a, 0)] can be decom- 
posed into different symmetry-adapted combinations of 
the basis functions, as shown in Table [j], where we have 
included, for the two-dimensional representations E g and 
E u , only one of the two possible partners. Note that we 
are thus able to essentially fit the eigenfunctions as well 
as the eigenvalues by this use of the block-diagonalization 
procedure. In the above listing, the first column is sim- 
ply the conventional band label, the second the conven- 
tional group symbol (g and u are the usual positive and 
negative parity labels), and the third column lists the 
basis expansion functions belonging to that particular ir- 
reducible representation. 

To begin the fitting procedure, it is necessary to make 
some choice of initial parameters. In our experience, the 
fitting procedure was not sensitive to the choice of ini- 
tial parameters, provided that the block-diagonalization 
procedure was used as the fit was optimized. 

For this reason, we have typically chosen a very sim- 
ple initial set of parameters. We usually set A = 1, e° 
near the expected center of the relevant band and the 
other e l a = for i = 1, 2, or 3, and the diagonal ele- 
ments of s —1. All of the off-diagonal elements for both 
the Hamiltonian and overlap matrices are started at zero. 
With this starting set, we typically optimize one of the 
cubic structures for 3 or 4 volumes spanning the range 
of volumes of interest, and allow the parameters to opti- 
mize with block- diagonalization. This provides the real 
starting TB parameters. 

The fitting process involves standard nonlinear least 
squares algorithmsEZl. Because this type of optimization 
process is prone to becoming trapped in local minima, the 
initial step described is important in order to provide a 
good starting point for the full fitting process. 

An alternative starting procedure has been described 
by NRLta that involves algebraically solving a simpli- 
fied nearest-neighbor tight-binding model at high sym- 
metry points for cubic crystal structures and then com- 
paring the solution to the LAPW bands at the specified 
points to determine a resulting reduced set of parame- 
ters. These are then used as the starting values for the 
nearest-neighbor-shell non-hybridizing matrix elements, 
namely {ssa, ppa, ppir, dda, ddn, dd5, ffa, ffn, ff8, ffcj>}. 
Other matrix elements are started with zero value. We 
have also tried this, more complicated, procedure. How- 
ever, we converged to fits of comparable quality as our 
standard approach (described above). 

In the full optimization we include fee, bec, and sc 
crystal structures in addition to the ground-state crystal 
structure (if different) and any other structures of inter- 
est. Because we know the correct block-diagonalization 



at high symmetry points and lines for the cubic phases, 
keeping these in the fit prevents the minimization from 
wandering into unphysical parts of parameter space. For 
all structures and volumes included in the fit, we use 
both the energy eigenvalues of the relevant bands as well 
as the total energy. During the minimization procedure, 
we typically weight the total energies by a factor of 1,000 
to 10,000 times that of individual energy eigenvalues, and 
include a penalty factor for singular overlap matrices. No 
block-diagonalization was done for the non-cubic struc- 
tures. 



B. The Database 



TABLE II: Data used in TB fits, including number of bands 
fit in each volume, Nb, and the number of volumes fit in each 
crystal structure, Nv- Also listed are the average root-mean- 
square errors in the fitted energy bands, (et), and the average 
error in the fitted total energy, (Et)- 



Element 


Structure 


N v 


N t 


<*6> [eV] 


(E T ) [eV] 


Be 


fee 


9 


3 


0.73 


0.0063 




bee 


7 


3 


0.69 


0.0053 




sc 


6 


3 


1.01 


0.0034 


Cu 


fee 


12 


6 


0.11 


0.0018 




bec 


8 


6 


0.12 


0.0021 




sc 


9 


6 


0.12 


0.0027 


Nb 


fee 


7 


6 


0.25 


0.0021 




bec 


8 


6 


0.25 


0.0018 




sc 


5 


6 


0.29 


0.0015 


U 


fee 


6 


9 


0.42 


0.0105 




bee 


5 


9 


0.36 


0.0075 




sc 


6 


9 


0.59 


0.0082 




oC4 


14 


18 


0.39 


0.0173 



To determine the TB parameters, we have fit to 
a database of highly accurate total energies and en- 
ergy bands from a series of full-potential linearized aug- 
mented plane-wave (FLAPW) calculations including lo- 
cal orbitalslia. The data for which each element was fit 
consisted of energy bands and total energies for a series 
of volumes in various crystal structures. Table [0] lists the 
number of volumes and crystal structure types for each 
of the elements studied here, Be, Cu, Nb, and U. 

Also shown in Table || is the number of energy bands 
fit for fc-points in the IBZ (47 IBZ /c-points were used for 
fee and bec lattices, 56 for sc, and 24 for a-U), which 
were also used in determining the total energies using 
a temperature broadening special point integration (us- 
ing a broadening of 2mRy). The block diagonalization 
procedure was carried out at all high-symmetry points 
as well as the midpoints of all high symmetry directions 
in the BZ, which made it necessary to break down the 
FLAPW energy bands by their irreducible representation 
at all such points. A representative plot illustrating the 
quality of the TB fit for low lying energy bands is shown 
in Figure [IB, in this case for fee U near the equilibrium 
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volume; note that we show bands for the foe structure 
rather than the ground-state structure in order to re- 
duce the complexity of the figure (a similar figure for 
a-U would have 2 times as many bands). The fit quality 
is quite good for a substantial energy region extending 
far above the fermi level. The highest energy bands in- 
volve higher-lying orbitals that are not included in our 
tight-binding basis, and hence the TB fit breaks down at 
high energy. For one crystal structure it is often possible 




terials (application to a threejjcenter sp 3 basis has been 
explored, for example, in Ref.Eil), and apply it to heavy 
elements. 



III. APPLICATION TO METALLIC ELEMENTS 

We now apply the TB total energy method to a range 
metallic elements from the periodic table. In order to 
provide a fair test for the TB method as an interpolation 
scheme, we have selected Be to represent light metals, 
Cu and Nb for transition metals, and U as a heavy f- 
electron metal. Table ID list the structural properties 
predicted by our TB fits, compared with FLAPW and 
experimental results. In general, the TB fits are very 
accurate, and highly transferable. The equilibrium vol- 
ume, i>o, and bulk modulus, B , as well as the pressure 



derivative of the bulk modulus, B' a , were determined by 
fitting the various equation p£ state curves to a second 
order Birch equation of stateEi. 

Beryllium. The alkaline-earth metal Be belongs to a 
part of the periodic table known for nearly free electron 
behavior, yet the properties of this metal, including the 
c/a ratio, are strongly influenced by the directional char- 
acter of the interatomic bonding. In determining our TB 
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FIG. 2: Electronic energy bands for fee U, a = 8.2, show- 
ing the quality of the TB fit (solid symbols) compared to 
the FLAPW results. Note that the first 9 bands were fitted 
throughout the IBZ. Low lying bands are fit more accurately 
since they also contribute to the total energy. Bands far above 
the fermi energy were not fit. 

to fit the energy bands so well that the difference between 
the LAPW and the TB bands are almost indistinguish- 
able. However, when all of the different structures are 
added into the TB fitting procedure, the overall quality 
of the fit to the individu al en ergy bands tends to degrade 
to that shown in Figure II B. In practice, we have found 
that improving the total energy fit lower than 1 mRy 
requires some compensating increase in the individual 
energy band errors. Ideally the two sets of errors should 
collectively decrease (after all, the total energy is the sum 
of the occupied band energies), but there may be some 
shortcoming in the functional form of the TB parameters 
that cause the two sets of errors to compensate for one 
another rather than for them to behave collectively. We 
hope to explore this issue in future work by considering 
alternate functional forms. 

These FLAPW results matched very well with a recent 
ab-initio study of the structural properties of U using 
FLAPW and a complementary technique using a basis 
set of Gaussian type orbitalsE3. We have neglected spin- 
orbit coupling for U, which has little effect on the bulk 
propertiesEl In a future paper we will extend this TB 
technique to include spin orbit coupling for f-electron ma- 
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FIG. 3: Equation of state curves for various crystal phases 
of Be, both in TB (solid symbols), and in FLAPW (solid 
curves) . 

fit for Be, we used a sp 3 d 5 basis with parameters fit to 
a database of 22 different volumes in the fee, bec, and 
sc phases. It was necessary to restrict the fit to only the 
three lowest energy valence bands (Ny in Table |J) f° r 
Be, as the bands corresponding to the atomic 3s states 
lie just above the fermi level, and our basis set has only 
a single s-type basis function (a fit using an additional s 
state would be able to fit higher energy bands). We found 
that an accurate model required inclusion of the d or- 
bitals, which are quite close energetically to the occupied 
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TABLE III: Structural parameters of hep Be, fee Cu, bee Nb, and a-U in the TB model , FLAPW calculation, and experiment 
(at a temperature where consistent structural and elastic data was available, noted in parentheses). The equilibrium volume is 
denoted by vo, and the bulk modulus and its pressure derivative by Bo and B' . y is the atomic position parameter in the oC4 
a-U structure 



method 


v (A ) 


7~> //nT)., \ 

.tfo(GPa) 


Bq 
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b/a 
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hep Be 














TB 


8.04 


127 


2.0 


1.565 






FT.APW 

J? 1 ii \ I vv 


7 QO 


123 


o.oo 


1 






lixpt(lOOK) 


8.05 


1 i ri) 

115 




1.571 






fee Cu 














1 D 


11.9 


i on n 
139. U 


0.1 








FLAPW 


11.9 


140.3 


5.1 








Expt(293K) c 


11.79 d 


138 


5.3 








bec Nb 














TB 


18.1 


160 


3.9 








FLAPW 


18.1 


161 


3.8 








Expt(120K) 


18.02 6 


171.& 










a-U 














TB 


20.24 


131 


2.8 


1.732 


2.044 


0.1017 


FLAPW 


20.14 


142 


5.0 


1.741 


2.073 


0.0990 


Expt 9 


20.58 


136 h 




1.734 


2.063 


0.1023 



a Rof. 22 

6 Ref. 23 
c Rcf. 2J 
d Ref. 25 
c Ref. 2E 
/Ref. 22 
9Ref. 2B 
ft Ref. 29 



s and p orbitals. The ground state structural parameters 
for the Be TB model are compared with experiment in 
Table [IF Figure || shows that our TB model successfully 
predicts the ground state crystal structure (hep) in addi- 
tion to an accurate prediction of the equilibrium volume 
and c/a ratio. The c/a ratio was determined by using a 
cubic polynomial fit to total energy calculations for 8 c/a 
values around the experimental value. 

Copper and Niobium. Equation of state plots compar- 
ing our TB model for the transition metals Cu and Nb 
are shown in Figures ^ and ||. The predicted structural 
parameters are given in Table For both transition 
metals the lowest six energy bands were fitted through- 
out the IBZ. For both elements the ground state structure 
is cubic, with a small energy difference between the face 
and body-centered structures, which the TB fit repro- 
duces quite well. Note the excellent agreement between 
the TB model and the FLAPW calculations, as well as 
the broad range of volumes and energies fitted. 

Uranium. For U, it has been shown that the bulk struc- 
tural properties are rather insensitive to ther-txeatment of 
spin-orbit coupling for the valence electronsEJ, so this el- 
ement makes a good test case for the scalar relativistic 
treatment that we have thus far implemented in our TB 
method. For U, the lowest 9 energy bands per atom were 
fit throughout the IBZ for the cubic crystal structures as 
well as the a-U structure (oC4) . The atomic position pa- 
rameter, y, was determined by performing a cubic fit to 




10 15 20 25 30 

Volume/Atom [A 3 ] 

FIG. 4: Equation of state curves for various crystal phases 
of Cu, both in TB (solid symbols), and in FLAPW (solid 



a six calculations of the total energy versus y. The ratios 
b/a and c/a were then found by using this value of y and 
performing a set of total energy calculations at the ex- 
perimental volume for twenty values of b/a and c/a. The 
ratios that minimized the total energy were then found 
using a 10 parameter, two dimensional cubic polynomial. 

The equation of state for the a, fee, bee, and sc phases 
are shown in Figure ^, which shows excellent agreement 
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Volume/Atom [A 3 ] 



FIG. 5: Equation of state curves for various crystal phases 
of Nb, both in TB (solid symbols), and in FLAPW (solid 
curves) . 

between the calculated total energies between FLAPW 
and TB, thus justifying TB as an accurate interpolation 
scheme, even for f-electron materials. Again the agree- 




a-U 



Q 

16 18 20 22 24 26 

Volume/Atom [A 3 ] 

FIG. 6: Equation of state curves for various crystal phases of 
U, both in TB (solid symbols), and in FLAPW (solid curves). 

ment between the TB and FLAPW calculations is excel- 
lent. 



IV. CONCLUSIONS 

We have presented the results of a study of the struc- 
tural properties of the metallic elements Be, Cu, Nb, and 
U using a tight-binding total energy method. Fits for the 



TB method were obtained from highly accurate FLAPW 
calculations. The TB models thus derived are highly ef- 
ficient, due to the treatment of valence electrons only, 
and highly accurate, proving to be transferable to new 
structures not included in the original fit. 

We have also shown that the NRL total energy tight- 
binding parameterization is highly accurate for f-electron 
as well as for d-electron systems. In addition we have 
demonstrated that forcing the tight-binding fitting pro- 
cedure to use the correct wave-function symmetry when 
fitting to specific energy eigenvalues at high symmetry 
k points and lines in the Brillouin zone leads to highly 
transferable tight-binding paramcterizations. This pa- 
rameterization seems insensitive to the initial guess used 
to start the least-squares minimization procedure. 

Transferability is a key issue: do the tight-binding cal- 
culations accurately reproduce what a high-quality first- 
principles method would calculate? Like the answer for 
any approximation, our conclusions are mixed. The re- 
sults of this paper show that the current NRL param- 
eterization make it possible to fit the total energy as a 
function of volume (for a significantly large range of vol- 
umes) for many different crystal structures. This is a 
very significant degree of transferability. 

With respect to other properties, additional detailed 
results will be published in future publications for appli- 
cations to specific materials. However, it is possible and 
relevant to make some pertinent observations here. Like 
any large parameterization, the properties that can be 
accurately calculated depend on how well one has sam- 
pled all aspects of the parameterizations in the data base 
used to fit the parameters. Because of the highly complex 
and nonlinear relationship between the parameters and 
the properties, it is always difficult to guess ahead of time 
how well one has sampled the full parameter space. If the 
relevant sampling has been done well for the properties 
of interest, those properties will be accurately predicted 
and compare well with first-principles calculations for the 
same properties. 

Some preliminary results indicate that elastic con- 
stants often are of about 10-20% accuracy. Transferabil- 
ity can be improved by specifically including in the fit 
finite distortions of the unit cell that sample those elastic 
constants (the distortions should be large enough to give 
a significant energy difference on the order of the differ- 
ence between other crystal structures life fee and bec, yet 
small enough so that the energy change is quadratic in 
energy with the size of the distortion so that one remains 
in the elastic regime) . Similarly, preliminary phonon cal- 
culations are of a similar degree of accuracy without any 
additional fine tuning of the parameters, but can be sig- 
nificantly improved by adding in a few frozen phonons of 
finite amplitude into the fit. 

In addition, it is worth pointing out that adding in 
the diamond lattice into the fits seems to improve the 
phonons all by itself (before adding in frozen phonons, 
for example). The diamond lattice seems to sample other 
parts of the parameterization that have not been sampled 
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by the other cubic crystal structures (e.g., with only four 
nearest- neighbors, it improves the sample of Eq. ( |13| ) to 
smaller p) . As more experience is gained, it may be possi- 
ble to find other crystal structures or other properties to 
include in the fits that similarly improve transferability. 
This is an ongoing work. 

Finally, the NRL parameterization is certainly not the 
only possible parameterization. Significant improvement 
of transferability may be achieved by finding better func- 
tional forms for the hopping, overlap, and diagonal ma- 
trix elements. 

For /-electron systems, because the rare-earth series 
of elements have localized / electrons (beyond Ce), the 
interesting systems are the actinide series. For this series 
it will be important to include spin-orbit coupling due 
to the large charge of the nucleus, especially for Np and 
Pu. We hope to explore addition of spin-orbit coupling 
to the tight-binding method for f-electron materials in 
future work. 

To conclude: based on the excellent results of the tight- 
binding total energies that we have found for a wide 



variety of single-element metals and crystal structures, 
the tight-binding approach appears promising for pro- 
viding highly accurate calculations of complicated ge- 
ometries such as defect states and distorted structures, 
where many atoms are required in large supercells, and 
first- principles methods are prohibitively expensive. 
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